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SUMMARY 

For steady multidimensional convection, the QUICK scheme has several 
attractive properties. However, for highly convective simulation of step pro- 
files, QUICK produces unphysical overshoots and a few oscillations, and this 
may cause serious problems in nonlinear flows. Fortunately, it is possible to 
modify the convective flux by writing the "normalized" convected control-volume 
face value as a function of the normalized adjacent upstream node value, devel- 
oping criteria for monotonic resolution without sacrificing formal accuracy. 
This results in a nonlinear functional relationship between the normalized var- 
iables, whereas standard methods are all linear in this sense. The resulting 
Simple High Accuracy Resolution Program (SHARP) can be applied to steady multi- 
dimensional flows containing thin shear or mixing layers, shock waves, and 
other frontal phenomena. This represents a significant advance in modeling 
highly convective flows of engineering and geophysical importance. SHARP is 
based on an explicit, conservative, control-volume flux formulation, equally 
applicable to one-, two-, or three-dimensional elliptic, parabolic, hyper- 
bolic, or mixed-flow regimes. Results are given for the bench-mark purely con- 
vective oblique-step test. The monotonic SHARP solutions are compared with 
the diffusive first-order results and the nonmonotonic predictions of second- 
and third-order upwinding. 


INTRODUCTION 

Successful modelling of strong convection is one of the most challenging 
problems in computational mechanics. If the truncation error terms in the 
numerical approximation contain second-order spatial derivatives (as in the 
case of first-order upwinding), simulated results are artificially diffusive 
and often grossly inaccurate. Central difference methods introduce propagat- 
ing numerical dispersion terms (odd-order derivatives) which may corrupt large 
regions of the flow with unphysical oscillations. Contrary to common belief, 
the extent of these oscillations actually increases for higher order (central) 
methods. Higher order upwind schemes have been successful in eliminating arti- 
ficial diffusion, while minimizing numerical dispersion. In the case of 
second-order upwinding (ref. 1), leading truncation error is a (potentially 
oscillatory) third-derivative term; however, the fourth-derivative numerical 
dissipation is large enough to dampen the dispersion to some extent. Third- 
order upwinding, exemplified in the steady-state, control-volume case by QUICK 
(Quadratic Upstream Interpolation for Convective Kinematics), has a leading 
fourth-derivative truncation error term which is dissipative, but higher order 
dispersion terms may still cause overshoots and a few oscillations when 
excited by (what should be) nearly discontinuous behaviour of the convected 
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variable (ref. 2). Currently many general-purpose elliptic solvers (replacing 
those previously based on variations of essentially first-order upwinding, 
such as older versions of the well-known TEACH code (ref. 3), for example) are 
now using either second-order upwinding (refs. 4 and 5) or QUICK (refs. 5 to 
20) as the basis for their convective transport solver. 

QUICK, in particular, has several attractive properties: no numerical 

diffusion (leading truncation error is fourth-order "dissipation" as distinct 
from second-order "diffusion"); low dispersion (the leading dispersion term is 
a small fifth derivative, strongly damped by the fourth-order dissipation); 
inherent convective stability (due to the upwinded curvature terms, even in 
the absence of physical diffusion); algorithmic simplicity (based on a conserv- 
ative control -volume flux formulation); and computational efficiency (in terms 
of total "cost" for a prescribed accuracy). QUICK also has excellent pressure 
prediction capability in Navi er-Stokes codes; in particular, computed stagna- 
tion pressure remains constant in isentropic regions (as it should), whereas 
this is not the case with other convection codes (ref. 6). Some groups using 
"QUICKened" TEACH codes together with TEACH's standard tridiagonal line solver 
have experienced convergence problems with QUICK in strongly recirculating 
flow simulations (ref. 21). However, this appears to be due to using a single 
sweep direction; with alternating-direction tridiagonal (or pentadiagonal ) 
line sweeps, QUICK is extremely robust and reliable under all flow conditions 
(ref. 22). Certainly, the explicit time-marching solution method described 
here presents no problems, even in the inviscid limit. 

QUICK'S single shortcoming is its tendency, under highly convective condi- 
tions, to produce overshoots and possibly some oscillations on each side of 
(what should be) steps in the dependent variable when convected at an angle 
oblique (or skew) to the grid. Even in one-dimensional flow, QUICK produces a 
few oscillations upstream of a sudden jump in the convected variable, under 
high convection conditions. By contrast, second-order upwinding does not have 
this defect in one dimension; but, as seen later, this method too produces 
strong overshoots in two-dimensional oblique-step simulations. In some appli- 
cations, small overshoots and a few oscillations may be tolerable— merely rep- 
resenting inaccurate resolution of the discontinuity. More likely, however, 
nonlinear processes such as steepening in shock waves or the local behaviour 
of a computed diffusion or viscosity coefficient will feed back and amplify the 
oscillatory error, and may lead to catastrophic divergence (ref. 13). Current 
practice with codes based on QUICK seems to be to revert to adding artificial 
diffusion in an ad hoc manner in order to suppress overshoots. For example, 
first-order upwinding might be used for k and e equations, while QUICK is 
used for momentum and scalar transport (refs. 20 and 23). The penalty for this 
"patch-up" procedure is not immediately obvious; but, given its poor track 
record, one should always be suspicious of first-order upwinding. 

Clearly a code retaining QUICK'S desirable attributes while eliminating 
unphysical overshoots and oscillations would be of great practical signifi- 
cance. Sharp monotonic resolution of thin shear layers, species density jumps, 
temperature discontinuities, shock waves, and other frontal phenomena is a fun- 
damental goal of computational fluid dynamics. The following sections of this 
paper will show how it is possible to modify the multi-dimensional QUICK scheme 
to achieve this goal while retaining QUICK'S third-order global accuracy and 
good stability characteri si tics ; and, perhaps surprisingly, this can be done 
with very little additional computational cost, because the standard QUICK 
algorithm (or a slight variation thereof) is used throughout the overwhelming 
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bulk of the flow domain — i.e., any (more expensive) modification is used only 
in thin regions requiring special treatment, thus representing only a small 
fraction of the overall number of grid points. 

The next section describes the Normalized Variable Diagram (NVD) — a plot 
of the locally normalized convected control-volume face variable with respect 
to the normalized adjacent upstream node variable. In this plane, standard 
methods such as first- and second-order upwinding, second-order central differ 
encing, and QUICK, are all represented by (different) straight lines. It will 
become clear that in order to satisfy both high accuracy and monotonicity, a 
nonlinear functional relationship is necessary. The choice of this nonlinear 
function is not unique. However, a simple scheme based on exponential upwind- 
ing has all the desired properties and is highly compatiable with QUICK (both 
use the same grid nodes for interpolation); hence this is used as the basis 
for the resulting Simple High Accuracy Resolution Program. Details of the 
development of exponential upwinding are given in the second section, where it 
is seen that this nonlinear scheme is also third-order accurate. A quantita- 
tive criterion is devised for deciding when to use the standard QUICK scheme 
(in "smooth" regions) and when to invoke the more sophisticated monotonic 
interpolation; this quite naturally depends on the normalized curvature of the 
convected variable. Because exponential upwinding does not cover the entire 
range of the NVD, it is necessary to devise ad hoc extensions to match with 
QUICK at the extreme ranges; this is achieved via simple piece-wise linear con 
structions, resulting in what has become known as the Exponential Upwinding or 
Linear Extrapolation Refinement (EULER). This section closes with a sketch of 
the EULER-QUICK algorithm in one dimension, for clarity. The following sec- 
tion outlines the two-dimensional algorithm in detail, and shows how it can 
easily be extended to three-dimensional steady flow. Finally, results are 
given for the well-known bench-mark two-dimenesional pure-convection oblique- 
step test— probably the most severe test for any convection scheme. This is 
used for a direct comparison between classical first-order upwinding, second- 
order upwinding, QUICK, and SHARP. As expected, first-order upwinding is 
extremely artificially diffusive. Rather surprisingly, second-order upwiding 
exhibits quite strong overshoots at some convection angles. QUICK gives 
steeper resolution of the jump region, but generates angle-dependent over- 
shoots and some oscillations. By dramatic contrast, SHARP retains the steep 
resolution of QUICK — but remains absolutely monotonic. Its overall character- 
istics seem virtually insensitive to flow-to-grid angle. 


NORMALIZED VARIABLE DIAGRAM 
Definition of Normalized Variables 

Consider the variation of a convected scarlar <j>(x,y,z) along a direction 
normal to a control -volume (CV) face, as shown in figure 1(a). For the CV 
face convecting velocity direction shown, QUICK involves the two adjacent node 
values (<t>o and 4>o together with that at the next upstream node (<j>y) in 
modeling the convected CV face value, 4>f . Note that the labelling of node 
values— downstream (D), central (C), and upstream (U)--depends on the velocity 
direction, as of course does the choice of node for <j>y. Figure Kb) shows 
the same information in terms of the locally normalized variable. 
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( 1 ) 


~ <t> - 4>y 

* ~ 

Note particularly that, in terms of normalized variables, Id = 1 and |y = 0. 

For example, for QUICK, on a uniform grid, the convected CV face varia- 
ble is (ref. 2) 

•frf = 2 (< **D + < * > C ) “ 8 ( ^D “ 2(<> C + V (2) 

so, in terms of normalized variables 

| f = \ (1 + l c ) " l (1 " 2 $ c + 0) <3) 

or more conveniently 
QUICK: 

| f = 0.75 + 0.75($ c - 0.5) (4) 

It should be clear that if <j>f is a function of <j>o, <t>c» and <|>u, then 
the normalized variable, |f, is only a function of lc (since |q = 1 and 
lu = 0). This is the basis of the Normalized Variable Diagram (NVD), which is 
a plot of the functional relationship between the normalized convected face 
value, |f, and the normalized adjacent upstrem node value, |q. 


Linear Schemes 

Equation (3) shows that, for QUICK, the normalized variable diagram is a 
straight line passing through (0.5, 0.75) with a slope of 0.75. Other well- 
known schemes also have linear characteristics. For example, first-order 
upwinding requires, using the present notation (which takes account of flow 
direction), zeroth-order upwind "interpolation" 

$f = 4*0 (5) 

or, in terms of normalized variables, simply 
First-order upwinding: 

If = Iq (6) 

Similarly, second-order central differencing being independent of SGN(u n ), is 
simply the linear interpolation 


4v = o 


v 


(7) 
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which becomes, in terms of normalized variables, 
Second-order central : 


$ f = 0.75 + 0. 5($ c - 0.5) (8) 

And second-order upwinding, given by linear, upwind-biased extrapolation 



C 


2 *U 


(9) 


can be written 
Second-order upwinding: 


- f $ c (10) 

The linear NVD characteristics, equations (4), (6), (8), and (10) are 
shown in figure 2(a). The corresponding normalized interpolations are shown 
in figure 2(b) for a specific value of (< 0.5). Note that three of the 
characteristics pass through the point (0.5, 0.75), labelled Q. First-order 
upwinding passes through the origin, 0, and the point P at (1,1); but it 
passses wel 1 below Q. 

Characteristics passing through Q can be written 

$ f = 0.75 + S($ c - 0.5) (11) 

where S represents the slope of the line. Using the original (unnormalized) 
variables, these can be written in terms of the (upstream-weighted) curvature. 


<t> f = 0 . 5<4> D + <t> c ) - CF(<{> D - 24> c + 4> u ) (12) 


where CF is the curvature 

factor. 



Clearly, S = 2CF + 0.5, 

and in 

specific cases 


QUICK: 

CF = 

1 q 3 

8 ’ b " 4 

(13) 

Second-order central : 

CF = 

0, S = \ 

(14) 

Second-order upwind: 

CF = 

2 . b - 2 

(15) 


In terms of normalized variables, equation (12) becomes 

= 0.5(1 + $ c ) - CF( 1 - 2$ c ) (16) 

Note that, in general, any (nonlinear) functional relationship between <j>f and 
4>C passing through Q can be written in the form of equation (12) or (16) 
provided CF is taken to be a function of <j>c, rather than a constant. Also, 
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by making a Taylor series expansion about the CV face locations, it is not 
difficult to show that any (in general, nonlinear) characteristic: (1) pass- 

ing through Q is necessary and sufficient for second-order accuracy and 
(2) passing through Q with a slope of 3/4 is necessary and sufficient for 
third-order accuracy. This means that any scheme based on a characteri s i ti c 
which can be written in the form of equation (16), with CF = CF($q), is at 
least second-order accurate, and if 

CF(0. 5) = 0.125 — S(0. 5) = 0.75 (17) 

the scheme is third-order accurate. This, of course, correlates with equa- 
tions (13) to (15), which show that the simple second-order schemes pass 
through Q with a slope other than 3/4, whereas the third-order QUICK scheme 
indeed has S = 3/4. Note that first-order upwinding cannot be written in the 
form of equation (16), since it does not pass through Q. The nonlinear NVD 
characteristic to be developed in the next section will pass through Q with 
a slope of 3/4, thus maintaining formal third-order accuracy. 

In the appendix, it shown that linear NVD characteristics which pass 
through the second quadrant may produce unphysical oscillations in steady one- 
dimensional convection. This is a well-known failing of central differencing, 
and may also occur to some extent with QUICK under high convection conditions 
(ref. 2). From figure 1(a), one sees Immediately that these two characteris- 
tics indeed pass through the second quadrant. Experience has shown that such 
schemes are also oscillatory in two-dimensional steady-flow simulations. 
Characteristics which pass through the fourth quadrant (i.e., below 0) are 
artificially diffusive. Thus, in order to avoid oscillations without being 
artificially diffusive, one necessary condi tion for the nonlinear characteris- 
tic to satisfy is that it must pass through the origin, 0. Numerical experi- 
mentation has also shown that NVD characteristics which passs above P are 
oscillatory in two-dimensions (although not necessarily so in one dimension— 
second-order upwinding being the classic example). Similarly, passing below 
P gives artificially diffusive results. So another necessary condition for 
the nonlinear characteri s i ti c is that it must pass through P. Behaviour of 
the nonlinear scheme to be developed can be summarized for the monotonic 
regime (0 < <j>c < 1 ) : 

The nonlinear NVD characteristic should pass through 0, P, and Q, with 
a slope of 3/4 at Q. 

For values less than 0 or greater than 1, the characteristic should 
be extended in a continuous manner, ultimately approaching the QUICK line for 
extreme values. The next section outlines the development of a scheme which 
satisfies the above criteria. 


EXPONENTIAL UPWINDING OR LINEAR EXTRAPOLATION REFINEMENT 
Exponential Upwinding 

Quadratic upstream interpolation is based on assumed local behaviour of 
the form 


<t> = a + bE + cE 2 


(18) 
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where % is a local spatial coordinate normal to the CV face, positive in 
the direction of the convecting velocity, as seen in figure 1(a). Evaluating 
a, b, and c, in terms of the local node values results in 


QUICK: 


4>(£) = <j>Q + 



i + 


% - 24> c ♦ » L 
V 2 Ax 2 


(19) 


and, of course, setting £ = Ax/2 results in equation (2) for 4>f. 

Now consider an entirely different type of interpolation through the same 
three node values, based on an assumed (upstream-weighted) exponential of the 
form 



4><0 

= A + B exp(CC) 

(20) 

Eval uati ng 

the three parameters in 

terms of node values leads to 



<j>(0)=<t>^ = A + B 

(21) 


4>(-Ax) 

= 4» u = A + Be‘ CAx 

(22) 


4>( Ax ) 

= 4> d = A + Be CAx 

(23) 

from which 

it is easily found that 

A = • 

(Vu “ *c) 

(24) 


% - 2 ^c + <t> u ) 


and 



4> f = A ± y(4> D - A ) ( 4>jj - A) 


(25) 


or, in terms of normalized variables 



v* c (i - v 3 - 1 
(1 - 2| c ) 


(26) 


with no ambiguity of sign on the square-root. This represents the desired 
Exponential Upwinding (EU) characteristic for the normalized variable diagram. 
Note that it is defined only in the monotonic regime (0 < < 1). There is 

an indeterminacy at = 0.5; but it is easy to show, using L'Hopital's rule, 
that 


$ f (0. 5) = 0.75 


(27) 
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(28) 


Similarly, straightforward differentiation results in 


3*: S( EU) * 0.75 

l/ eu 



0.5 


showing that Exponential Upwinding is tangent to the QUICK line at = 0-5. 
This means, of course, that Exponential Upwinding is third-order accurate. 

The exponential upwinding characteristic shown in figure 3(a) in relation to 
QUICK (dashed), while figure 3(b) shows the corresponding normalized EU and 
QUICK interpolations for a small positive value of with the correspond- 
ing <j>f values shown in figure 3(a). 


Modification Criterion 

Note that the EU curve lies quite close to QUICK over a fairly wide 
range, near ~ 0.5. This suggests a very simple modification strategy for 
deciding whether to use the basic QUICK scheme or the more sophisticated 
interpolation: 

If | 0 . 5 - $^1 < const, use QUICK (29) 

where, from figure 3(a) a value of const = 0.15 might be considered reasona- 
ble. Multiplying by 2, this becomes 

If | 1 - 2$ c | < 0.3, use QUICK (30) 

or more directly, when written in terms of unnormalized variables 

If 1^ - 2<t> c + <|> D | < 0.3|<j> D - <^1 , use QUICK (31) 

Since the left-hand side of this equation is proportional to the curva- 
ture (normal to the CV face) of the convected variable, the modification cri- 
terion is a quantitative statement of the desire to use QUICK in "smooth" 

(i.e., small-curvature) regions of the flow domain. This will be the case in 
the bulk of the flow, since high curvature (rapid change in gradient) occurs 
only in thin regions involving a small number of grid points. Thus, although 
Exponential Upwinding is more expensive than QUICK, it is only used in a small 
fraction of the computational domain (if at all), so that the overall strategy 
is extremely cost-effective. 


Nonmonotonic Regime 

Since Exponential Upwinding is only available in the monotonic regime, 
the question remains as to the best procedure to adopt for <j>c > 1 and 
<t>C i 0- Eor <j>Q above 1, a simple and apparently robust strategy is to use 
the continuous extension 

<j>f = $ c for 1 < <i> c < 1 .5 (32) 

returning (again without loss of continuity) to QUICK at <j>c > 1-5. Note that, 
although this portion of the overall nonlinear NVD characteristic happens to 
coincide with that of first-order upwinding, it does not degrade the order of 
the overall algorithm— which is determined solely by equations (27) and (28). 
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The whole concept_of "order" based on Taylor series only has meaning for smooth 
behaviour, near $0 -> 0.5 (i.e, vanishing curvature). Statements such as those 
often made in relation to TVD schemes (ref. 24) (that such schemes are only 
first-order accurate near extrema: $ 0 = 0 , or 1 , in the present notation) are 

totally meaningless and can be quite misleading. What is actually meant is 
that, in the present notation, $f = $q for $c < 0 and $c _> 1 , which happens 
to coincide with the first-order upwinding characteristic in the nonmonotonic 
regime, but not (necessarily) near $q -»• 0.5. 

The negative-$c regime requires somewhat more care in designing an exten- 
sion from Exponential Upwinding, which ends at (0,0). A characteristic which 
rejoins QUICK at some finite negative $ 0 -value seems desirable; this could be 
done with a straight-line characteristic through ( 0 , 0 ) with a slope less than 
3/4. However, it i sjmportant to avoid a certain critical point on the QUICK 
characteristic (at $c = - / 3/2), since, as shown in the appendix, traversing 
this point on the QUICK line could lead to unphysical oscillations under cer- 
tain circumstances. It is thus better to rejoin the QUICK line below the crit- 
ical point. This is adequately accomplished by the ad hoc straight-line 
characteristic 

$f = | $C for -1 <. $0 1 0 (33) 

continuing along the QUICK line for $0 < -1 . The complete composite NVD char- 
acteristic is shown in figure 4, representing an Exponential Upwinding or 
Linear Extrapolation Refinement of QUICK. 


The EULER-QUICK Algorithm 

The one-dimensional algorithm is summarized here for reference for each 
CV face: 


(i) 

Designate upstream and downstream nodes on the basis of 

SGN(u n ) 

(ii) 

If |$ D - $y| < 10 -5 (say), use QUICK, otherwise 


(i i i ) 

Check if equation (31) is satisfied (this will account 
bulk of the flow field). 

for the 

( i v) 

If not, compute $ c = (4> c - ^y) / (<t> D - <^), and find $ f 

by 

(V) 

QUICK, equation (4), if $^ < -1 , or $^ > 1.5, or 
0.35 < $£ < 0.65, 


(vi) 

$^. = 0.375 $£ if -1 < $q 0, 


(vii) 

Exponential Upwinding, equation (26) for 0 < $^ < 0.35, 
0.65 < $ c < 1 , or 

and for 

(vi i i ) 

$^ = $ c for 1 < $ c <1 .5. 


( i x) 

Then, reconstruct the (unnormalized) face value, 
<t>f = 4>U + (4 >q _ 4>y)$f 




In this way, the convective flux at the left face can be computed 

CFLUXL(i) = CXL( i ) *<j> f (34) 

where CXL(i) = u(i) = u ( i ) At /Ax is the local Courant number at the left face 
for CV(i). Because of the conservative CV formulation, the right-face flux 
at CV(i) is just the left-face flux at CV(1+1), regardless of velocity direc- 
tion. The explicit update algorithm for pure convection then becomes, quite 
simply 

<{>" +1 = <t>" + CFLUXL(i) - CFLUXL( i + 1 ) (34) 

In the steady state, of course, 4> l ? +1 = (j) 1 ?, so that the FLUX terms must 

balance. Diffusive flux terms are treated in an analogous fashion (using 
DXL = Tj, At/Ax 2 , r being a diffusion coefficient or viscosity): 

DFLUXL(i) = DXL(i)-(4> i - 4*^) (36) 

involving the simple linear difference across the left face, which is consist- 
ent with the third-order treatment of the convective fluxes (ref. 2). Control- 
volume averaged source terms can be added, if appropriate. 


MULTI-DIMENSIONAL ALGORITHM: SHARP 

Two-Dimensional QUICK Scheme 

Figure 5 shows a two-dimensional control volume, with attention focussed 
on convection across the left face. In the situation shown, uq, is positive 
to the right, so for a control volume centered at node (i,j), the following 
designation of upstream and downstream variables results, in the direction nor- 
mal to the face 

4>q = <t>(i ,j) ; 4> c = <t>(i-1 , j) ; 4> u = 4>(i-2,j) for u ft > 0 (37) 

and the upstream-weighted transverse nodes are 

<j> T = $0-1 ,j+l); <j> B = 4>( i - 1 , j-1 ) for u^ > 0 (38) 

It should be clear which nodes would be involved for uq, < 0. 

For the basic two-dimensional QUICK scheme (ref. 25), the convected value 
averaged at the left face is (for ug, > 0) 

♦d = 2 (<t> D + V ■ 8 (<t *D “ 2 *C + V + 24 (<i> T " 2<j) C + < * > B ) <39) 

which includes the linear interpolation term, the upstream normal curvature 
term, and a small term representing the effect of (upstream-biased) transverse 
curvature, in computing the face average. Note that the first two terms are 
identical with the one-dimensional formula, equation (92). Thus, one possible 
strategy for the two-dimensional SHARP is to follow the one-dimensional algo- 
rithm, steps (i) to (ix), above, and then simply add the appropriate transverse 
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curvature term (depending in the sign of uj,) . One needs to store left-face 
(convective plus diffusive) fluxes 

FLUXL(i.j) = CXL( 1 ,j ) *<J>^ - DXL(i ,j)*[<j>(1 ,j) - 4>( 1-1 , j ) ] (40) 

and bottom-face fluxes 

FLUXB(i.j) = CYB(i ,j).<(> b - DYB(1 JM*(1 ,j> - 4>< i , j-1 > 3 (41) 

The explicit update algorithm is then extremely simple, using an over-writing 
assignment statement. 

Set: $( i ,j ) = <t>(i,j) + AtS*(i,j) + FLUXL( i ,j ) - FLUXL(i+l,j) 

+ FLUXB(i.j) - FLUXB( i , j + 1 ) (42) 

where S* is the CV-averaged source term. Note that, in terms of storage 
requirements, it is necessary to store two flux arrays for each transport 
variable. 


Curvature Based Algorithm 

An alternative strategy, that can be used for one-, two-, or three- 
dimensional simulations, is to write, for the average left-face convected 
value (suppressing j and k indexes, for convenience) 

<j> a = \ (<t> i + <{>._.,) - CF-CURVN + (CURVTY + CURVTZ) (43) 

where CURVTY and CURVTZ are the respective upstream-weighted transverse 
curvature terms in the other coordinate directions, and the upstream-biased 
normal curvature can be written 

CURVN = ^ + i + 4*1-2 - SGN(u^) *(4>^ + ^ - 3<|>. + 3<J>. 

(44) 


which takes account of velocity direction automatically. For the standard 
QUICK algorithm, of course, the normal curvature factor, CF, is a constant 
(= 1/8). The explicit SHARP algorithm can be implemented simply by writing, 
from equation (16). 


$ f ($ r ) - 0.5(1 + $ r ) 

CF = CF($ r ) = — (45) 

L (2$ c - 1) 

which is shown in figure 6, along with the nonmonotonic extensions, correspond- 
ing to the complete NVD characteristic of figure 4. 

Again, the algorithm begins by computing |<$>d - <t>u|, and going immediately 
to QUICK if this is less than a specified small number. If not, one computes 
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in the usual way, and then goes to algebraic expressions representing the 
behaviour shown in figure 6 for CF($q), as follows: 


For $£ < -1 or <j>£ 2 1*5: CF = 0.125 

(0.5 + 0.125 $ r > 

For -1 < $ c < 0: CF = — 

(1 - 2$ ) 

C 

U$ c - 1) 

For 1 < $ r < 1 • 5 : CF = - — 

L (2$ c - 1) 

For 0.5 < < 0.7, use the following quadratic approximation: 

CF = 0.125 - 0.2609($ c - 0.5) + 0.13613($ c - 0.5) 2 


(46) 

(47) 

(48) 

(49) 


And finally, for 0 < <j>c < 0.3, and 0.7 < 4 >q < 1» use the exact Exponen- 
tial Upwinding formula, combining equations (26) and (45). 


$ r - (1 + $ r )($ r - 0.5) - V$ r (l - <j> r )3 

CF = — w—— ~ — (50) 

(1 - 2 $ c )^ 

It should be clear that most of the grid points will be in the smooth region 
(<j>C - 0.5) and thus will involve equation (49). This empirical formula is 
graphically indistinguishable from the exact EU curve in this region. 

A less expensive strategy uses the following approximate formulas, in 


addition to equation (46): 

For -1 < $ c < 0: CF = ^ + ($ c - 1) (51) 

For 0 < $ c < j : CF = \ - | J$~ c (52) 

For 1 < $ c < 1 : CF = 1 (1 - $ c > (53) 

For 1 < $ c < | : CF . 1 - i (| - » c ) <54) 


In this case, most grid points will involve the simple linear formula given in 
equation (53). Note that CF = 1/8 when <j>Q = 1/2, thus maintaining third- 
order accuracy even in this simple approximate scheme. The main difference 
between equations (47) to (50) and equations (51) to (54) is that in the latter 
case the monotonic extensions are approximated by straight lines for CF($q), 
rather than segments of hyperbolas; since this is an ad hoc procedure in either 
case, the exact shape of these extensions is immaterial, provided they revert 
to the QUICK value (CF = 1/8) in a reasonable manner (and avoid the critical 
point) . 
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Extension to Three Dimensions 

The three-dimensional algorithm parallels the two-dimensional version, 
simply adding the upstream-biased transverse curvature for the third direc- 
tion, and then setting up (convective plus diffusive) fluxes for "left" (L), 
"bottom" (B), and "far" (F) faces of each CV cell. The explicit update then 
becomes 

Set: 

4><i,j,k) = <|>( i , j ,k) + AtS*( i , j ,k) + FLUXL(i,j,k) - FLUXLC i + 1 , j ,k> 

+ FLUXB(i,j,k) - FLUXB( i+1 , j ,k) + FLUXF(i,j,k) - FLUXF( i , j ,k+l ) (55) 

where, as usual, flux consistency is guaranteed by the control -volume form- 
ulation; e.g., the "near" face flux, FLUXN( i , j ,k) , has been replaced by 
FLUXF ( i , j , k+ 1 ) . In the three-dimensional case, of course, it is necessary to 
store three flux arrays for each transport variable. 


Numerical Boundary Conditions 

As a general principle, QUICK boundary-condition treatment is used 
(ref. 25); i.e., quadratic extrapolation normal to the boundary to set up 
external pseudo-node values, using the given physical conditions together with 
enough interior node values to perform the extrapolation. This works well for 
all transport variables unless there is a discontinuity oblique to the grid 
very close to the boundary. In such a case, the classic solution of local mesh 
refinement can be used. Alternatively, other high-resolution forms of local 
behavior can be developed on a case-by-case basis. For example, boundary- 
layer flow with a no-slip condition and strong pressure gradient can be mod- 
eled on a relatively coarse grid in a manner compatible with QUICK (and SHARP) 
by assuming a local behaviour normal to the wall such as 

u = C^ 17 " + C 2 y (56) 

i.e., a power-law-plus-linear profile, with an assumed value of n. This is 
really a three-point interpolation, since the condition u(0) = 0 has already 
been assumed; it would replace the usual (three-point) quadratic variation 
(used in QUICK) 

u = c 1 y + c 2 y 2 (57) 

again assuming u(0) = 0. In each case, the coefficients are computed using 
two internal node values normal to the wall. Equation (56) allows much more 
rapid variation in the CV cell adjacent to the wall; it can be used in the 
same spirit as logarithmic wall-functions (ref. 3). Naturally, details regard- 
ing convective and diffusive fluxes need to be worked out in individual cases; 
but this is a straightforward matter. 

At inflow boundaries it is sometimes convenient to set up two external 
pseudo-nodes so that the usual (internal) QUICK or SHARP algorithm can be used 
directly. Figure 7 shows a situation in which 4 >bc is given at an inflow 
boundary. Local quadratic behaviour in terms of the distance from the bound- 
ary implies 
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4>4 " 4 > 2 ^ f$A ~ 2< ^ + ♦ 


4 T y 2 ,2 




( 58 ) 


2 Ax 


thus 


^BC = 4" 2 " ^3 4 (<i> 4 “ < ^ > 2 ) + 8 (<<> 4 ^3 + < * > 2 ) 


(59) 


Solving for 4 ^ gives 


8 1 
*2 = 3 ^BC " 2<l> 3 + 3 ^4 


(60) 


And since for a quadratic function, the third-difference is zero, <!>, is given 
by 


4>1 = 34>2 _ 34>2 + $4 


(61) 


For outflow conditions, one may assume zero-local-curvature in the flow 
direction, unless other physical conditions are specified. This is clearly 
not as restrictive as zero-local-gradient, and may allow the use of smaller 
computational domains for a given accuracy. 


Time-step Restrictions 

If one makes a classical von Neumann analysis of the (one-dimensional) 
QUICK scheme, assuming unsteady convection and diffusion in an infinite 
domain, the resulting time-step restrictions are rather stringent (ref. 2). 

In fact, the convection-controlled restriction is formally the same as the 
simple forward-time-central-space requirement (ref. 27) on the Courant number, 
c = uAt/Ax, 


c < i- (62) 

A 

where = uAx/r is the grid Peclet number (or Reynolds number). This would 
imply, of course, that purely convective flow (P^ - ®) could not be simulated 
by the explicit QUICK scheme. Fortunately, the formal von Neumann analysis 
does not apply to steady-state algorithms on a finite grid. In particular, 
since the instability indicated by violating equation (62) is at the long- 
wavelength end of the Fourier spectrum, the imposition of a long-wavelength 
cutoff (corresponding to a finite grid) results in a much less restrictive 
condition (ref. 27). This can be written in an accurate simplified form as 
(ref. 25) 

? 2 

c < 5 “ + (63) 

r A 2N^ 

where NAx = \*, the cut-off wavelength, from which it can be seen that even 
in the "inviscid" limit (P^ - ®) there is always a nonzero time step available 
for explicit solution of the QUICK algorithm on a finite grid. Even this can 
be violated when using time-marching toward a steady-state solution, because 
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any unstable modes tend to be suppressed by the steady-state boundary condi- 
tions. Finally, the nonlinearity of the SHARP scheme allows further violation 
of formal time-step restrictions. As a rule of thumb, one should try Courant 
numbers about one order of magnitude larger than that suggested by equa- 
tion (63). One other point needs to be noted, however: when simulating 

purely convective flow, short-wavelength transients may "rattle around" in the 
computational domain. And although these are not unstable, there is no mecha- 
nism for their damping, and a general level of "background noise" is thereby 
generated. In real flows, of course, there will always be some physical damp- 
ing, although this might be extremely small. Thus, when simulating "purely 
convective" flow, the following strategy is suggested: include modeled diffu- 

sion terms (finite P^) , but increase P^ until the solution is independent 
of its value. A nominal value of P^ = 10^ satisfies this requirement, while 
completely quieting transient noise. This magnitude is certainly far beyond 
any value likely to be encountered in practice; i.e., it represents a computa- 
tional "infinity." 


Variable Grids 

For clarity, the development of SHARP has been based on the assumption of 
a uniform grid in each coordinate direction. Various levels of generalization 
are possible. For example, in two dimensions it is a simple matter to extend 
the formulas to a uniform rectangular grid, Ax = const, Ay = const * Ax; and 
similarly in three dimensions. This is merely reflected in the definition of 
individual component Courant numbers (and diffusion parameters). The next 
level of generalization involves locally expanding (or contracting) rectangu- 
lar grids, with "expansion ratios" such as r x = Ax ^ + ] / Ax ^ , etc. In princi- 
ple, one could set up analogous formulas for 4>f incorporating r x , ry (and 
r z in three dimensions), as has been done for QUICK (refs. 2 and 25) using a 
different notation. However, it turns out that, for QUICK, simply using the 
constant-grid-spacing formulas on a variable grid results in negligible errors, 
provided the adjacent mesh-width ratios lie within the range 0.8 ~ 1.25, i.e., 
up to approximately a 125 percent local expansion ratio (ref. 4). This gives 
a wide range of flexibility in designing variable rectangular meshes without 
going to the added complexity of variable-grid formulas. Since this conclu- 
sion was based on a Taylor series expansion about control-volume faces, and is 
therefore related to QUICK'S NVD behaviour near = 0.5, the same conclusion 
must be reached regarding SHARP. The extension to nonrectangul ar quadrilat- 
eral grids is not yet documented; however, it is reasonable to speculate that 
"mild" distortion would not result in significant errors, even though the for- 
mal order of accuracy is reduced. 


TEST PROBLEM RESULTS 
The Obi ique-Step Test 

Figure 8 shows the well-known bench-mark test problem consisting of pure 
convection of an upstream transverse step profile in a scalar field imposed at 
the inflow boundaries of a (square) computational domain, in this case, 

Ax = Ay = ^ (64) 
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There are two additional rows of pseudo-nodes upstream of the inflow bounda- 
ries, and one additional set downstream. The convecting velocity is of 
constant magnitude and flows at the same angle 9, oblique to the grid, every- 
where. The location of the boundary step is chosen so that the exact conven- 
ed step passes through the mid-point of the grid, for reference. Note that <j> 

= 0.5 along the step itself, whereas, <j> = 1 everywhere above and <j> = 0 below, 
as indicated. 

For 0 = 45°, the exact solution is shown in orthographic projection in 
figure 9. Note that this particular computer-plot routine interpolates line- 
arly between specified grid-point values. Figure 10 gives the results for 9 
= 45° using classical first-order upwinding for all CV face fluxes. Clearly, 
this is grossly in error due to the artificial cross-grid diffusion inherent 
in this method. A quantitative indication of the error is given by 

ERROR = ^ | ^computed - ^exact I (65) 

summed over all computed (interior) grid points; the magnitude is noted in the 
figure captions in each case. Figure 11 gives the corresponding results for 
second-order upwinding. Although the main rise is considerably steeper than 
first-order, and the ERROR much smaller, the most obvious feature is the (anti- 
symmetrical) pattern of overshoots — of considerable magnitude! This is a seri- 
ous problem, and one that needs to be addressed by research groups propounding 
the use of second-order upwinding as a general-purpose convection solver 
(refs. 4 and 5). In fact, second-order upwinding, at this angle, gives larger 
overshoots than QUICK, seen in figure 12. The QUICK and second-order-upwind 
results are qualitatively similar; but note that QUICK'S step resolution is 
considerably steeper, and compare the quantitative ERROR magnitudes. Finally, 
for this 45° angle, figure 13 gives the SHARP results — essentially the same 
steep resolution as QUICK, but with the overshoots "clipped off" (and smoothed) 
to give absolutely monotonic simulation. Note the drop in ERROR relative to 
QUICK. 

Of course, 45° is a special angle (actually, the worst case for first- 
order upwinding in terms of ERROR magnitude), so two other angles will be 
considered 


9 = artan 



34° 


( 66 ) 


and 

9 = artan » 56° 

The results— again for first-order upwinding, second-order upwinding, QUICK, 
and SHARP— are shown in the final series of figures, as indicated. In each 
case, first-order upwinding is artificially diffusive, second-order upwinding 
and QUICK are oscillatory with much lower ERROR, but SHARP always gives uni- 
formly steep and monotonic results, essentially independent of flow-to-grid 
angle . 


DI SUSS ION AND FORECAST 

SHARP represents a new generation of multidimensional monotonic convective 
solvers of high formal accuracy (in this case, third order). Other similar 
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schemes can be constructed by devising alternate nonlinear characteristics in 
the normalized variable diagram. For example. Gas ke 1 1 and Lau's Sharp Mono- 
tonic Algorithm for Realistic Transport (SMART) is based on a piece-wise lin- 
ear characteristic (ref. 28) consisting of the QUICK line for the bulk of the 
monotonic range in <j>Q, but deviating via ad hoc straight-line segments to 
pass through (0,0) and (1,1) in the NVD. SMART and SHARP give virtually iden- 
tical results for critical steady two-dimensional pure-convection problems 
such as the oblique-step test. These schemes are similar in some respects to 
certain types of flux-limiter and so-called "TVD" schemes developed for simu- 
lating shock phenomena in inviscid compressible flows. In fact, there is a 
one-to-one correspondence between the normalized variable diagram and the 
flux-limiter versus gradient-ratio diagram discussed, for example, by Sweby 
(ref. 29). It is not difficult to show that the flux-limiter factor, FLF, 
(Sweby's cp), is related to the normalized variables used here by 


FLF 




and that the gradient ratio, r, is given by 


(68) 


r 


<1 - * c > 


(69) 


Note that the important region near <j>Q = 1 is banished to large (positive or 
negative) r values in Sweby's diagram, whereas the behaviour of the 4>f(<j>c) 
characteristic is immediately obvious in the NVD used here. In terms of unnor- 
malized variables. Equation (68) is simply 


FLF 



V 

V 


(70) 


where the numerator is the difference between the modelled face value, 4>f, and 
first-order upwinding, <j>Q, and the denominator is the difference between 
second-order central differencing, <}>fCEN = i/2(4>q + 4>q) , and first-order 
upwinding. In unsteady flows, one can obtain the central-difference time- 
averaged (Lax-Wendroff ) CV face value from 





»c> 


(71) 


where c is the Courant number. If one makes the same (second-order time- 
accurate) assumption for the nonlinear face value, averaged over time, 

<<|> f > = 4> f - c(<|> f - <|> c ) (72) 

then equation (70) can be written, cancelling the factor (1 - c) 


FLF 


<< 4 > f > - 4 > c > 

- 4 > c > 


(73) 


17 


which, of course, is the basic definition of the flux-limiter factor. This 
means that many of the flux-limiter schemes previously developed for unsteady 
(one-dimensional) gasdynamics will be applicable to steady multidimensional 
highly convective flows, as well, using the conservative control -vol ume formu- 
lation described here. 

It now appears possible to obtain even better resolution of discontinui- 
ties by using "nonlinear" schemes of higher-order accuracy, in the sense of 
using more than three grid points in estimating local behaviour (normal to 
control-volume faces). Whether this is based on flux-limited higher-order 
polynomial schemes or more sophisticated forms of (nonpolynomial) interpola- 
tion, the strategy of using a robust scheme, such as third-order upwinding, in 
the bulk of the flow domain and switching to the (presumably) more costly com- 
pution only where necessary (in thin layers) will remain highly cost-effective. 
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APPENDIX 


Consider a linear characteristic in the NVD of the form 


^ = + i 


(A. 1 > 


In one-dimensional simulations, if oscillations occur, they will take the form 
of an alternating geometric series, decaying upstream: <j>i = 1 , _ 1 = -K, 

<t>i_2 = etc. Let 4 >d = 1 ; then <j>£ = -K, <f>(j = K^. The corresponding nor- 
malized variable is 


r* *U _ -K - K 2 _ -K 

C 1 - K 2 1 ' K 


(A. 2) 


where the "star" signifies a critical value of <j>c (corresponding to the oscil- 
latory behaviour). The accompanying face value is <j>f* = 0, as seen in 
figure A. 1(a); thus the normalized face value becomes 



This gives a quadratic equation for K, the appropriate root being 

I, _ S - Vs 2 - 41(1 - S - I~ 

2(1 - S - I) 

For example, for QUICK (S = 0.75, I = 0.375) 

K(QUICK) = 2 V3 - 3 = 0.4641 
and the critical point is given by 

$ C *(QUICK) = - V3/2 = -0.8660 


(A. 3) 


(A. 4) 


(A. 5) 


(A. 6) 


and 


$*(QUICK) = - | <V3 - 1) = -0.2745 (A. 7) 

as shown in figure A. 1(b). 

Note that for second-order central differencing (S = 0.5, I = 0.5), 

K = 1, implying undamped oscillation. Also, if the characteristic passes 
through 0 (I = 0), then K = 0; i.e., perfect step resolution. Finally, if 
I < 0 (i.e., the characteristic passes through the fourth quadrant), K is neg- 
ative, implying a nonosci 1 latory geometric decay, or artificial diffusion. It 
should be clear that any linear characteristic passing through the second quad- 
rant (I > 0) will have a corresponding critical point in the third quadrant. 
More general results concerning nonlinear characteristics can be determined. 

In particular, any (in general, nonlinear) characteristic passing through the 
second quadrant will have a potentially oscillatory critical point in the 
third quadrant. This can be seen by imagining a local first-order Taylor 
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expansion of the form of equation (A.l) about candidate critical points; pro- 
vided the trajectory enters the third quadrant from the second, there will 
always be values of S and I satisfying equations (A. 2) to (A. 4) at some 
point on the trajectory on the third quadrant. Trajectories entering the 
third quadrant through 0 or the fourth quadrant may also have critical 
points, as would the scheme shown dashed in figure A. Kb), which rejoins the 
QUICK scheme above its critical point. Clearly, the most satisfactory design 
is to pass through 0, staying low enough to avoid potential critical points. 
This is the case with the ad hoc linear extension chosen in constructing the 
EULER-QUICK scheme. 
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(A) ORIGINAL VARIABLES. <B) NORMALIZED VARIABLES. 

FIGURE 1. - NODE VARIABLES IN THE VICINITY OF A CV FACE (AT DASHED LINE). 
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(A) NORMALIZED VARIABLE DIAGRAM. (B) NORMALIZED INTERPOLATIONS 

(FOR O c < 0.5). 

FIGURE 2. - FIRST-ORDER UPWINDING, SECOND-ORDER UPWINDING, SECOND-ORDER CENTRAL 
DIFFERENCING, AND THIRD-ORDER UPWINDING (QUICK). 



(A) NORMALIZED VARIABLE DIAGRAM. (B) NORMALIZED INTERPOLATIONS. 


FIGURE 3. - EXPONENTIAL UPWINDING (SOLID CURVES) SHOWN IN RELATION TO QUICK (DASHED). 
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FIGURE 12. - OBIIOUE-STEP TEST RESULTS FOR THIRD-ORDER UPWINDING (QUICK). 0 » R5°. ERROR - 16.6. 
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SECOND-ORDER UPWINDING, 6 = 56° 


FIGURE 21. - OBLIQUE-STEP TEST RESULTS FOR SECOND-ORDER UPWINDING, 0 = 56°. 


ERROR * 28.7. 



FI6URE 22. - OBLIQUE-STEP TEST RESULTS FOR 


OUICK, 6 s 58°. ERROR = 23.4. 
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servative, control-volume flux formulation, equally applicable to one-, two-, or 
three-dimensional elliptic, parabolic, hyperbolic, or mixed-flow regimes. 

Results are given for the bench-mark purely convective oblique-step test. The 
monotonic SHARP solutions are compared with the diffusive first-order results 
and the nonmonotonic predictions of second- and third-order upwinding. 
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